# References:
# [1] Some Product Formulas for Theta Functions in One and Two Variables, 
#     by David Grant
# [2] Theta Functions and Singular Torsion of Elliptic Curves,
#     by John Boxall and David Grant
# [3] Resultants of Division Polynomials II: Generalized Jacobi's Derivative 
#     Formula and Singular Torsion on Elliptic Curves,
#     by David Grant
from ThetaFunction import thetaDeriv, theta, thetaSecondDeriv
from EtaFunction import eta
from numpy import pi
from EllipticCurve import EllipticCurve
from sympy import Integer as Int
from sympy import div
from ThetaProds import derivProd, prod
from sympy.abc import y
from SigmaFunction import sigma

# Ref [1], Equation (1), page 223
def test1(m, e):
    print prod(m, e.tau)
    print (-1)**(m-1) * m * eta(e.tau)**(m**2-1)

# Ref [1], Lemma 2 (i), page 229
def test2(e):
    print thetaDeriv(0.5,0.5,0.0,e.tau)
    print -pi * theta(0.0,0.0,0.0,e.tau) * theta(0.5,0.0,0.0,e.tau) * \
        theta(0.0,0.5,0.0,e.tau)
    print -2 * pi * eta(e.tau)**3

# Ref [1], Theorem 1 (i) and (ii), page 233
def test3(e):
    for (a,b) in [(0,0),(0,0.5),(0.5,0)]:
        print a,b
        print theta(a,b,0.0,e.tau)**8 * derivProd(3,e.tau,a,b)
        print ((-1)**(2*a+1) * 2**8 * pi**8/(3.0**5)) * eta(e.tau)**32
    print derivProd(4, e.tau, 0.0, 0.0, [(2,0),(0,2),(0,0),(2,2)])
    print (-(pi**12)/8.0) * eta(e.tau)**36

def test4a(n,e,z):
    f = e.divPoly(n)
    print sigma(n*z, e.w1, e.w2)/(sigma(z, e.w1, e.w2)**(n**2))
    print f.eval(e.pFunc(z))

# Ref [3]
def test4(n,e,z):
    f = e.divPoly(n)
    c = -2 * pi * eta(e.tau)**3
    print c**(n**2-1) * theta(0.5,0.5,n*z,e.tau)/\
        (theta(0.5,0.5,z,e.tau)**(n**2))
    print f.eval(e.pFunc(z))

# Ref [3]
def test5(n,e,a,b):
    f = e.divPoly(n)
    g = e.secondDeriv(f)
    c = -2 * pi * eta(e.tau)**3
    z = (a*e.tau+b)/float(n)
    fac1 = n**2 * c**(n**2-1)
    tz = theta(0.5,0.5,z,e.tau)
    tpz = thetaDeriv(0.5,0.5,z,e.tau)
    tpnz = thetaDeriv(0.5,0.5,n*z,e.tau)
    tppnz = thetaSecondDeriv(0.5,0.5,n*z,e.tau)
    print g.eval(e.pFunc(z))
    print fac1 * (- 2*n*tpz*tpnz + tppnz*tz)/(tz**(n**2+1))

# Ref [3]
def test6(n,e,a,b):
    f = e.divPoly(n)
    g = e.secondDeriv(f)
    z = (a * e.tau + b)/float(n)
    c = -2 * pi * eta(e.tau)**3
    print 2 * c**(n**2) * n**3 * (-1)**(b+1) * \
        thetaDeriv(0.5 + float(a)/n, 0.5 + float(b)/n, 0.0, e.tau)/\
        (theta(0.5 + float(a)/n, 0.5 + float(b)/n, 0.0, e.tau)**(n**2+1))
    print complex(g.eval(e.pFunc(z)))

# Ref [3]
def test7(n,e):
    f = e.divPoly(n)
    g = e.secondDeriv(f)
    print ((-1)**((n-1)/2) * Int(2)**(n**4-1) * \
           (Int(1)*pi)**(n**4-n**2)*(Int(1)*eta(e.tau))**\
           (2*n**4-3*n**2+1)*n**(2*n**2-4)*derivProd(n,e.tau)).evalf(10)
    print (f.resultant(g)**2)/(n**(n**2+1))

def test8(n,m,e):
    fn = e.divPoly(n)
    fm = e.divPoly(m)
    if m % 2 == 0:
        (fm, r) = div(fm, 2*y)
    if n % 2 == 0:
        (fn, r) = div(fn, 2*y)
    print fn.resultant(fm)
    print e.discrim()

e = EllipticCurve((0.5+1.2j),1,True)
#e = EllipticCurve(20.0,30.0)
#e = EllipticCurve(1.3,2.2)
#z = (1.25+0.j)
#test4(7, e, z)
test5(5, e, 6, 7)
